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ABSTRACT 

We study numerically the axisymmetric relativistic Bondi-Hoyle accretion of a su- 
personic ideal gas onto a fixed Schwarzschild background spacetime described with 
horizon penetrating coordinates. We verify that a nearly stationary shock cone forms 
and that the properties of the shock cone are consistent with previous results in New- 
tonian gravity and former relativistic studies. The fact that the evolution of the gas is 
tracked on a spatial domain that contains a portion of the inner part of the black hole 
avoids the need to impose boundary conditions on a time-like boundary as done in 
the past. Thus, our approach contributes to the solution to the Bondi-Hoyle accretion 
problem at the length scale of the accretor in the sense that the gas is genuinely en- 
tering the black hole. As an astrophysical application, we study for a set of particular 
physical parameters, the spectrum of the shock cone vibrations and their potential 
association with QPO sources. 
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1 INTRODUCTION 

Bondi-Hoyle accretion is the process of evolution of an 
infinite uniform wind of gas moving near a central object, 
or conversely a compact object moving on a uniform 
gas with constant velocity ( |Bondi fc Hoyle 1944[ ). This 
process is associated to various phenomena related to gas 
processes near stars or compact objects. It shows interesting 
properties when considered within the Newtonian and rela- 
tivistic regimes, which have been explored based on several 
numerical studies in the classical regime ruled by Newto- 
nian gravity (|Shima et al. 1985] IMastsuda et al. 19871 
Fryxell fc Taam 1988| ISawada et al. 19891 

Mastsuda et al. 19911 IMastsuda et al. 19921 

Ruff'crt fc Ar nett 19941 IRuffert 1994al IRuffert 1994bl 
Rufi'crt 19961 IRufi'ert 19971 IBenensohnT997l 

Nagae et al. 2004| [Blondin fc Raymer 20"l2| ), in which 
the most important subject are the consequences on the 
morphology of the wind and the supersonic shocks that 
develop; a summary of results under Newtonian gravity can 
be found in ( |Foglizzo et al. 2005[ ). 

Unlike in the Newtonian regime, the relativistic ap- 
proach allows the study of the Bondi-Hoyle accretion in 
regions where the gravitational field is strong, for in- 
stance, near the event horizon of a black hole. Some 
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studies have been pushed forward in this direction. The 
first one was performed by (|Petrich et al. 1989|) : in this 
work they carry out axially symmetric numerical simu- 
lations in order to study the different accretion patterns 
developed by the gas during the accretion onto a black 
hole. Later on, (|Font fc Ib ancz 1 998bllFont fc Ibanez 1998al 
IFont et al. 19981 IFont et~a l. 1999) revisited the results using 
high resolution shock capturing schemes that have shown to 
be much more accurate methods than those used in the past 
consisting in the addition of artificial viscosity. It is worth 
pointing out that in all these papers the morphology of sub- 
sonic and supersonics fiows was the most important sub- 
ject, because when the wind has supersonic velocity a high 
density shock cone forms on the opposite side of the wind 
source, and such cone has interesting properties like resonant 
oscillations or flip-fiop motion. On the other hand, in the as- 
trophysical context, a recent work by (|Donmez et al. 201 1|) 
dedicates important efforts to study a potential relation be- 
tween QPOs and the behavior of the gas density in the rel- 
ativistic Bondi-Hoyle accretion. 

Considering full general relativity, Bondi-Hoyle accre- 
tion has been studied by (|Farris et al. 20l'0|) in the con- 
text of the merger of supermassive black hole binaries. 
Also axisymmetric magnetohydrodynamics Bondi-Hoyle ac- 
cretion was studied in (|Penner et al. 2010[) . On the other 
hand (|Zanotti et al. 20111) studied the case of the relativistic 
Bondi-Holye accretion considering radiative processes and 
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the case of ultrarelativistic axisymmetric Bondi-Hoyle ac- 
cretion was recently considered in (|Penner et al. 2012|) . 

In this paper we present a relativistic study of super- 
sonic axisymmetric Bondi-Hoyle accretion onto a spheri- 
cally symmetric black hole described with horizon pene- 
trating coordinates. We assume the condition that the gas 
is a test fluid and does not distort the geometry of the 
space-time. Previous studies on the subject involve the use 
of coordinates that are singular at the event horizon, and 
therefore the evolution problem has to be set on a domain 
that requires an artificial inner boundary outside the black 
hole's horizon, which requires the implementation of efflcient 
boundary conditions there, specially when such boundary is 
close to the event horizon where the state variables tend 
to diverge, in order to prevent errors to propagate toward 
the numerical domain and consequently affect the numerical 
calculations. What the use penetrating coordinates offers, 
is the possibility to define the inner boundary inside the 
black hole's horizon and furthermore to avoid the imple- 
mentation of boundary conditions there, because the light 
cones using penetrating coordinates remain open and ma- 
terial particles move naturally toward inside the black hole. 
In other words, the fluid is allowed to fall into the black 
hole. Some basic models of spherical accretion of dark mat- 
ter by supermassive black holes already use this condition 
HGuzman fc Lora 2011allGuzman fc Lora 2011b[) . 

Even though our approach represents an improvement 
to the complete study of Bondi-Hoyle accretion, it is only 
an alternative to solve the problem at the inner boundary 
related to the length scale of the accretor. At a different 
scale, the accretion radius is defined to approximately decide 
when a particle falls into the compact object or not, and 
such scale is determined by the velocity of the wind and the 
equation of state of the gas. Traditionally, the Bondi-Hoyle 
problem is attended in this second scale, by assuming the 
accretor is point-like, whereas in our case we are exactly 
doing the opposite. Numerically it seems that one has to 
choose between these two regimes, that is, it is possible to 
apply wind accretion to astrophysical processes if realistic 
wind velocities are considered, and consequently the accretor 
length scale is small enough as to be considered a point- 
like source that cannot be numerically resolved; conversely, 
if one wants to resolve the accretor (like in our case) the 
accretion radius length scale is unresolved unless the velocity 
of the wind is assumed to be high. This difficulty of studying 
both scales at the same time is the main reason why the 
relativistic wind accretion has not been fully resolved at the 
moment. 

The association of numerical results on Bondi-Hoyle ac- 
cretion onto black holes to astrophysical observations in- 
volves a series of parameters that in principle should be as- 
trophysically motivated. For instance the adiabatic index 
and internal energy of the gas, and the velocity of the wind 
if one wants to infer the mass and angular momentum of the 
black hole, or viceversa; furthermore, properties of the gas 
Eissociated to heat transfer of cooling processes expected to 
happen would involve an extra set of parameters that should 
be observationally motivated. Perhaps the most solid bounds 
among all these parameters is the relative velocity between 
the accretor and the wind, which is of the order of 100-1000 
km/s in binaries and may achieve from a few thousands of 
km/s ([Gonzalez et al. 2007|) in superkicked black holes re- 



sulting from the collision of two black holes with coUimated 
emission of gravitational radiation up to very recent results 
indicating 15000 km/s ( |Sperhake et al. 2011[ l in similar sce- 
narios. 

The paper is organized as follows. In section [2] we de- 
fine the relativistic Euler equations describing the fluid onto 
a fixed Schwarzschild black hole in Eddington-Finkelstein 
coordinates, whereas in section [3] we explain the numerical 
methods we use to solve the relativitic Euler equations on 
a curved background space-time; in section [4] we show the 
morphology and mass accretion rates for the accretion of 
an axisymmetric homogenous wind and in[S]we analyze the 
frequencies associated to the rest mass density oscillations 
inside the shock cone. Finally in [6] we present our conclu- 
sions. 



2 RELATIVISTIC HYDRODYNAMIC 
EQUATIONS 

For a generic space-time, relativistic Euler equations can be 
derived from the local conservation of the rest mass and the 
stress-energy tensor T^j^ of the fluid model considered: 

V,(T''") =0, 

V,{pou) = 0, (1) 

where po is the rest mass density, it'' is the four-velocity of 
the fluid and is the covariant derivative consistent with 
the four-metric of the space-time (jMisner et al. 1973|l . 
Notice that our calculations assume geometric units where 
G = c = 1. We consider the gas accreted by the black hole 
is a perfect fluid, which can be described by the following 
stress-energy tensor 

T^. =/9o/ii*''M"+TO^", (2) 

where p is the pressure and h the relativistic specific en- 
thalpy given by 

/i = l + e+— . (3) 

po 

where e is the rest frame specific internal energy density of 
the fluid. 

In order to track the evolution of the fluid it is con- 
venient to write down the relativistic Euler equations as 
flux balance laws (jMarti et al. 19911 [Banyuls et al. 1997| 
IFont et al. 2000|l . for which we require the space-time met- 
ric to be written in the 3-1-1 splitting approach of general 
relativity. In this paper we consider the accretion onto a 
spherically symmetric black hole for which the most general 
line element reads 

ds^ = —{a^ — 'yrrP^ P^)dt'^ + 2'yrrl3^ drdt + ^rrdr^ 

+ -1gg{de^ +smed4>''), (4) 

where a is the lapse function, the shift vector has only one 
component /3* = (/3'',0,0), 7ij are the components of the 
spatial metric associated with the space-like hypersurfaces 
and = {t, r, 6, cj}) are coordinates of the space-time that 
we choose to be spherical in the spatial part. 
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A fingerprint of our paper is tliat we use penetrat- 
ing coordinates, because tlie previous researcii uses mostly 
Scliwarzscliild coordinates, wliicli are singular at the event 
horizon, and most authors implement an artificial time-like 
boundary near the event horizon in a region where the grav- 
itational field is extremely strong and numerical bound- 
ary conditions in such region might not be as appropriate 
as in the far region. We instead use horizon penetrating 
Eddington-Finkelstein slices to describe the black hole, for 
which Q, /?' and read 



part of the 4- velocity u' , as = ^ + 77, where W is the 
Lorentz factor W = , 

The system of equations (f ^lOp is a set of four equa- 
tions for either the primitive variables po,v^ ,e,p or the 
conservative variables D, Jr, Je,T and p or e. Therefore it is 
necessary to close the system of equations. As usual we close 
the system using an equation of state relating e = e{po,p). 
In this work, we choose the gas to obey an ideal gas equation 
of state given by 



1+ 



2M 
r 



I _|_ ML I ' ^' ^ 



J- n , 2 2.2, 

7ij = aiag[l -\ ,r ,r sm ( 



(5) 



Since we consider the test field approximation to be valid 
in our numerical experiments, there is no need to evolve the 
geometry of the space-time, instead we keep our background 
space-time fixed and describe the evolution of the fluid on 
top of such space-time. The general equations obtained from 
IT} for a fluid with axisymmetry in spherical coordinates for 
metric (Q read 



h 



^/7 

(6) 

Here U is a vector of conservative variables, F* are the fluxes 
in the spatial directions and S is a vector of sources. These 
quantities are defined in terms of both the primitive vari- 
ables (po, v'~,v^ , e) and the conservative variables themselves 
as follows: 
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pohW'^ve 
pohW^ -p~ poW 



(7) 
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(8) 



(9) 



(10) 



In these expressions, 7 = det(^ij) is the determinant of the 
spatial 3-metric, V are the Christoffel symbols of the 
space-time and = (w'",u*,0) is the 3- velocity measured 
by an Eulerian observer and defined in terms of the spatial 



(r-i)po' 



(11) 



where F is the ratio of specific heats or adiabatic index. 
The relativistic speed of sound plays an important role in 
the definition of the initial conditions of the wind of gas we 
evolve, and for the equation of state can be written as 

^2 ^ pr(r-i) 

where its asymptotic value or its maximum permitted value 
is c. = y/V^. 



3 INITIAL DATA AND NUMERICAL 
METHODS 

3.1 Evolution 



We evolve the system of equations (|6ll0p us- 
ing a high resolution shock capturing method 
( [LeVeque 1992[). We specially use the HLLE Riemann 
solver (jHarten et al. 1983 fc Einfeldt 1988)l in combination 
with the minmod linear piecewise reconstructor. For the 
integration in time we used a second order Runge-Kutta 
integrator. 

We evolve the gas on the domain ([reaic, rmax] x [0, it]) x 
[0, 2n) in spherical coordinates using a two dimensional code, 
and a uniformly spaced grid along coordinates r and 9. 
The radial domain deserves a careful description. On the 
one hand, rexc defines a spherical boundary that we can 
choose inside the black hole's event horizon called an exci- 
sion boundary originally implemented to study the evolu- 
tion of black holes (jSeidel fc Suen 1992|> and recently incor- 
porated to the study of hydrodynamics (jHawke et al. 2005|l . 
That is, a chunk of the domain inside the black holes horizon 
is removed from the numerical domain from < r ^ Vexc in 
order to avoid the black hole singularity and the steep gra- 
dients of metric functions near there; since the light cones 
inside the horizon (located at r = 2M) point towards the 
singularity, there is no need to impose boundary conditions 
at r = rexc, instead, the fluid simply gets off the domain 
through such boundary. We choose the excision radius to be 
rexc = 1.5M, which is both far enough from the singularity 
at r = and provides a buffer zone 1.5M < r < 2M for the 
fluid inside the horizon to flow smoothly from the horizon to- 
wards the excision boundary. This is a rather improvement 
compared to previous results where Schwarzschild coordi- 
nates are used and the excision boundary is chosen outside 
the black hole, in a time-like region that requires the imple- 
mentation of boundary conditions. 

On the other hand r^nax defines an exterior spherical 
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boundary. We distinguish between two hemispheres called 
downstream and upstream as in (|Mastsuda et al. 19871 
IFont fc Ibafiez 1998b[l . The upstream hemisphere is defined 
as the part of the sphere where the gas is entering the do- 
main and the downstream is the hemisphere through which 
the gas leaves the domain. In the downstream hemisphere 
we impose outflow boundary conditions and in the upstream 
boundary we always consider the initial asymptotic value of 
all variables as set initially. 

Another important issue is that rmax has to be bigger 
than the accretion radius race (see below) if one wants shocks 
to form, and experience tells us that using rmax ~ lOracc 
provides adequate time and space scales to see the shock 
cone until a nearly stationary stage. 

At the axis defined by S = 0, tt the sources are extrap- 
olated from the cell next to the boundary and demand 
to be odd there. We also implemented an atmosphere, that 
is, we define a minimum value of po that avoids the diver- 
gence of the specific enthalpy and errors that propagate to 
the other variables. We set the density of the atmosphere to 
10"", which we found to allow accuracy and consistency of 
our numerical results. Even though in our numerical exper- 
iments we do not measure such a small rest mass density 
during the stationary regime, the implementation of the at- 
mosphere works only to prevent the density to be small in 
rarefaction zones that may form before the stationary stage. 
Finally, since the fiuxes depend both on the primitive and 
the conservative variables (see ((S]) and ((51) we reconstruct 
the pressure in terms of the conservative and the other prim- 
itive variables using a Newton- Raphson algorithm. 

We performed tests of our code for different radial do- 
mains and resolutions and found optimal number of cells 
depending on the different models, being particularly inter- 
ested in maintaining a sufficiently large domain preserving 
the spatial resolution within a convergence regime. Also in 
the appendix we present a set of more numerical standard 
tests for our numerical implementation. 



3.2 Initial data 

We set the wind to move along the direction of z initially, 
with constant density and pressure. We characterize the ini- 
tial velocity field in terms of the asymptotic initial velocity 
Hoc as done in ([Font fc Ibafiez 1998b|) : 



zVoo cos C^, 



zUoo sin 9, 



(13) 
(14) 



where the relation v — vtv^ — v^o is satisfied. With this 
conditions, the gas initially with constant rest mass density 
moves along the ^-direction filling the numerical domain. 

The initial density and pressure profiles are chosen in 
such a way that relation (|12|l is satisfied. In order to do 
this, we fix the initial density to a constant pini and give 
an asymptotic value for the sound speed , then using 
equation (|12|) the initial pressure can be obtained as follows 



A subtle point here is that we choose a value for to 
construct the initial pressure, with the only condition that 
in order to avoid negative or zero pressures we need the 
condition < y/V — 1 to hold. Finally, using the equation 
of state ITT} we calculate the initial internal energy ei„i. 

At this point Voa and Csoo are two important param- 
eters to set up the initial velocity field, and we find use- 
ful to define the relativistic Mach number at infinity in or- 
der to parametrize our initial data = Wvao /WaCsoc = 
W Mod/Ws, where W is the Lorentz factor of the gas Ws is 
the Lorentz factor calculated with the speed of sound and 
A^oo is the asymptotic Newtonian Mach number we use to 
parametrize our initial configurations. When this number is 
bigger than one it is said that the flow is supersonic an oth- 
erwise subsonic. In the supersonic case a high density shock 
cone forms and because the density in the cone oscillates the 
matter transport of such oscillations could be associated to 
QPO sources. 

A very important scale is the accretion radius defined 

by 



M 



(16) 



'pini — 



«(r-i) 



(r(r-i)-cLr)- 



(15) 



otherwise all the gas in the domain would be accreted and 
such domain would then be inadequate to observe a shock 
cone. This radius is the Bondi accretion radius, where M is 
the accretor mass, in our case a black hole of mass M = 1 
(|Petrich et al. 1989|) . 

The parameter space to study the Bondi-Hoyle accre- 
tion is enormous and the free parameters are the initial rest 
mass density of the wind pini, the value of F, 

^oo , Csoo or 

equivalently the internal energy eini of the gas or JVl . We de- 
cide to fix the initial rest mass density pi„i — 10^^ in units 
where the mass of the black hole is M = 1, which works 
in the test fluid approach we are assuming, although there 
are other precedents where a much higher density is used 
(|Donmez et al. 20111) . We also use this density because for 
a range between 10~* < pini < 10~* we find the shock cone 
and its oscillatory behavior. In Table [T] we present the pa- 
rameter space we explore. Another possible parameter would 
be the mass of the black hole, nevertheless we choose units in 
which Af = 1, and automatically spatial and time measures 
will be in units of M; then our results are mass independent, 
and in order to specialize to a particular astrophysical case 
where M is given in units of solar masses, we only need to 
rescale the units of space and time appropriately. We stress 
that we use parameters that allow the numerical track of 
the accretion process, which is restricted mainly by the size 
of the numerical domain; explicitly, if one assumes that ei- 
ther the gas or the black hole moves with speed of the order 
of Vaa ~> Csoo ~ lOOfcm/s, then Tacc/M ~ 10^, then as- 
suming the outer boundary of the numerical domain is at 
rmax ~ lOracc it would be ~ 10^ times the size of the black 
hole radius. In the most optimistic case, when a black hole is 
moving with a velocity of the order of 10^ fcm/s as a result of 
a superkick as response to the coUimated emission of gravita- 
tional radiation ( [Sperhake et al. 2011[ |, the accretion radius 
would be of the order of race ~ lOOM and rmax ~ lOOOM. 
This is at the moment a restriction to carry out accurate nu- 
merical calculations, because on the one hand one needs a 
considerably big domain and on the other one needs enough 
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Table 1. Set of parameters we use for our numerical experiments. 
We use different values of and the same density pi„i = 10~^ 
in all cases. The model in bold face was not solved numerically 
due to the computer limitations described in the text, instead it 
is an extrapolation with the other similar cases to Mach 1. We 
use the two values of F to test the shock formation and accretion, 
and use only F = 4/3 for the QPO model. 



Model 


F 




Moo 




Cs^ = 0.1 


Mia 


4/3 


0.5 


5 


3.84615 


Mu 


5/3 


0.5 


5 


3.84615 


M2a 


4/3 


0.4 


4 


5.88235 


Mib 


5/3 


0.4 


4 


5.88235 


M3a 


4/3 


0.3 


3 


10 


M3i 


5/3 


0.3 


3 


10 


Mia 


4/3 


0.2 


2 


20 


Ma 


5/3 


0.2 


2 


20 


Ms 


4/3 


0.1 


1 


50 


= 0.08 


Me 


4/3 


0.4 


5 


6.00926 


My 


4/3 


0.32 


4 


9.19118 


Ms 


4/3 


0.24 


3 


15.625 


Mg 


4/3 


0.16 


2 


31.25 


Mio 


4/3 


0.08 


1 


78.125 


Cs^ = 0.05 


Mil 


4/3 


0.25 


5 


15.3846 


Mi2 


4/3 


0.2 


4 


23.5294 


Mi3 


4/3 


0.15 


3 


40 


Mi4 


4/3 


0.1 


2 


80 


Ml5 


4/3 


0.05 


1 


200 



resolution as to resolve the black hole. What we do here is 
that in order to cover a numerical domain using the nec- 
essary resolution to achieve convergence, we require the do- 
main to be smaller than for the aforementioned astrophysical 
scenarios. We do this by considering rather high velocities of 
the gas that imply smaller values of race ~ lOM (see (|16|l 'l. 

3.3 Diagnostics 

In order to diagnose the physical quantities of the system 
in our simulations, we implement detectors located at vari- 
ous radii in the numerical domain, that is, we define spheres 
where we calculate scalar quantities. In particular, we track 
the mass accretion rate as a function of time. In order to 
do so we calculate the mass accretion rate on a sphere spe- 



cialized to the case of a spherically symmetric space-time in 
spherical coordinates 

Mace = -2n J Dr'^ ^v'' - sin6'd6', (17) 

at various spherical surfaces, including the black hole event 
horizon. The mass accretion rate helps diagnosing when the 
accretion has stabilized. 

The most important quantity we measure is the rest 
mass density of the gas along the axis inside the shock cone, 
that is at S = and at the location of the different detectors. 
This scalar is used to measure how the density oscillates in- 
side the shock cone. At post-step we perform a Fourier trans- 
form of this scalar in order to study the predominant oscil- 
lation modes of the density and study the relation between 
such oscillation and its potential relation to QPO sources. 



4 PROPERTIES OF THE SHOCK CONE AND 
MASS ACCRETION RATES 

Depending on the properties of the gas flow, a shock cone 
appears when the wind is supersonic. This shock cone is 
a region where the density is significantly higher than the 
density of the wind itself and it forms behind the black hole 
in the opposite side of the source of the wind. 

In Fig.m we show the morphology for the different mod- 
els we have considered in Table [1] As the parameter space is 
very large, we only chose two values for the adiabatic index 
r = 4/3,5/3. Moreover, we only consider four initial wind 
supersonic velocities of the gas at infinity corresponding to 
2, 3, 4, 5 Mach. We then show that the use of penetrating co- 
ordinates does not change the well known morphology found 
in previous relativistic studies when using time-like internal 
boundaries, that is, the bigger the Mach number the smaller 
the shock cone angle, and also when F increases the shock 
cone angle increases. 

In Fig. [2] the mass accretion rate for different values of 
the Mach number and F is presented. When the Mach num- 
ber and the adiabatic index increase, the system reaches the 
stationary regime faster. This has been already discussed in 
great detail by (|Font fc Ibafiez 1998b[) and to us represents a 
test. On the other hand, some morphological effects of using 
horizon penetrating coordinates for non axial accretion and 
flip-flop behavior of the shock have been recently presented 
(|Cruz-Osorio et al. 201211 . 

Once we have described some of the properties of the 
shock cone we measure the density at various points as a 
function of time. What happens is that the density vibrates 
as shown in Fig. (3] In our approach we deflne the inner 
boundary inside the black hole where the light cones are 
oriented toward the black hole singularity and are open as 
opposed to the BL coordinates that shows light cones that 
are closed when approaching the event horizon. Thus, the 
fluid particles move inside the black hole naturally. 

The fact that the shock cone vibrations are not asso- 
ciated to numerical artifacts, allows us to explore, as a toy 
astrophysical model, a possible relation between such oscil- 
lations inside the shock cone and QPOs frequencies. 
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Figure 1. Logarithm of the density that illustrates the morphology of the shock cones in various scenarios when the system has achieved 
a nearly stationary regime. In the top panel we present the cases for the wind in models (from left to right) Mia, M^a M2a and A/ia, 
whereas in the bottom for models M^f,, M^i,, M2b and A/ift described in Tabic [l] 
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Figure 2. We show the mass accretion rate Mace with respect to 
the proper time for different values of the Mach number and F. As 
we expect the accretion stabilizes faster when the Mach number 
and r increase. This quantity is calculated in a detector located 
at r = 2.1M, close to the event horizon. Despite the fact that wc 
are calculating the mass accretion rate near the event horizon, no 
high oscillations show up and our measurements are accurate. 



5 A POSSIBLE MODEL OF QPO SOURCES 

A rather appealing potential application of the vibrations of 
the density of the gas within the shock cone is related to the 



QPO type of emission. In (|Donmez et al. 201ip it is found 
the first work studying the possibility that the vibrations 
of the density within the shock cone could be associated to 
QPO frequencies. Some differences in our approach are that 
in our analysis uses penetrating coordinates in the descrip- 
tion of the black hole, whereas in (|Donmez et al. 201ip they 
are using singular coordinates and our study is axisymmet- 
ric whereas in (|Donmez et al. 201111 slab symmetry is used. 
With these tools we ask whether the most powerful peaks 
are approximately within the range of frequencies observed 
for QPOs in terms of the black hole mass. 

It is worth to mention that in our relativistic approach, 
aside of the restrictions of the realistic velocity of the wind, 
we do not include non-adiabatic cooling and heating pro- 
cesses since this is a first step toward an exhaustive and 
realistic analysis in the general relativistic regime truly al- 
lowing the gas entering the black hole. The first steps on 
that direction involving radiation terms have been recently 
published (|Zanotti et al. 20111 [Roedig et al. 2012[ ). 

We now analyze the oscillations as measured at vari- 
ous distances from the black hole. In Fig. [3] we show vari- 
ous aspects of the oscillations. Firstly we show the oscilla- 
tions as measured by a detector located at distance 12. IM 
in proper time; it can be seen that after an initial tran- 
sient the rest mass density approaches a nearly stationary 
regime; secondly we show a time window showing the sig- 
nal as measured by detectors located at various distances, 
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Figure 3. At the top we show how the density oscillates in time 
as measured by the detector at 12. IM. In the second plot we 
show a zoomed time window showing the signals as measured by 
various detectors equally spaced being the first one located at 
2.1 and the farthest one at 16.1; it can be seen the correlation 
among the signals measured by successive detectors; it is also 
shown that high frequency modes that are measured by detectors 
close to the hole are not seen by farther detectors. This result is 
complemented with power spectra below. 




500 1000 1500 2000 2500 3000 3500 4000 4500 5000 
t 

Figure 4. Self-convergence of the L2 norm of the density along 
the positive 2-axis, where the shock cone forms. The ratio among 
successive resolutions is 1.5, and the convergence factor is defined 
by l.S'?. The accuracy of our algorithms is expected to be first 
order due to the development of shocks. 



Figure 5. Spectrum of density oscillations in proper time for the 
model Mil, measured at diff'erent detectors located at: (Top- 
left) r = 2.1M, (top-right) r = 6.1M, (bottom-left) r = 12.1M 
and (bottom-right) r = IG.lAf . 



ranging from 2.1M to IG.lAf; from this plot we learn that 
detectors closer to the black hole measure oscillations dom- 
inated by high frequency modes, whereas detectors far from 
the black hole measure signals dominated by lower frequency 
oscillations; moreover, this plot shows that signals among 
detectors are well correlated and that we are not measuring 
only numerical noise. We support our calculations with the 
a self-convergence test of the density in Fig. 3] 

In order to interpret the oscillations of the density we 
perform a Fourier transform of the signals measured by dif- 
ferent detectors like those in Fig. [3] using the proper time at 
the location of each detector. The results appear in Fig. \^ 
where the high frequency modes are shown to dominate near 
the black hole and do not appear far from the black hole. 
We concentrate in modes that may be global and detected 
in a big part of the domain as suggested by the analysis 
in (|Donmez et al. 2011[l . As a particular case, we point out 
with an arrow in Fig. [S] a frequency peak that appears in 
all the detectors. That it appears in all our detectors is an 
indication that the corresponding mode is global and this is 
the reason why we choose it to correspond to the frequency 
in our analysis. In all the cases in our parameter space the 
spectra show the presence of such mode. This is the reason 
why in all the cases studied we consider the spectrum as 
measured with a detector located at r = 12. IM where the 
spectrum is very clean and the most powerful peak corre- 
sponds precisely to this mode. 

We then proceed to explore the chances that such fre- 
quency lies within observed QPO frequency ranges. Our pa- 
rameter space covers the models described in the Table for 
various values of the speed of sound and velocity of the wind. 
For each case we locate the most powerful peak as measured 
by the detector at r = 12. lA/ and obtain a proper frequency 
of oscillation of the gas within the shock cone. In this paper 
we push the wind velocity values down to 0.05c. We con- 
sider that the smallest velocity a wind requires to form a 
shock cone is that of the speed of sound. However, when 
the speed of sound is small, the accretion radius is big in 
terms of the length scale of the black hole (see the bold face 
parameters in the Table and eq. (|16|l ) and this imposes a 
computational limitation on the set of parameters that can 
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Figure 6. We show the frequencies of oscillation of the rest mass 
density in the shock cone in terms of the wind velocity for different 
values of the speed of sound. 



be used to complete a simulation, because one requires not 
only a big spatial domain, but also enough resolution near 
the black hole and in the whole domain. 

In Fig. [S] we show the frequencies obtained for the sim- 
ulations in our Table. Circles correspond to the results ob- 
tained from simulations considering a sound speed Csoo ~ 
0.1. Analogously, squares and triangles indicate our results 
with Csoo = 0.08 and Csoo ~ 0.05 respectively, being the 
model Mi5 an extrapolation to Mach 1 given the numerical 
domain restrictions described. 

In Fig. [7] we compare the frequency window we find in 
all our simulations with observational results. We ask about 
the value of our parameters needed for our frequencies to 
lie within the range of those observed from QPO sources. 
Each point in each of the diagonal lines defines a proper 
frequency of the peak we have selected in physical units for 
several different scenarios, in terms of the black hole mass. 
Each one corresponds to a different value of the velocity of 
the wind in units of the speed of sound. We obtain these 
lines by simply scaling our results, which were calculated in 
units of M = 1. In the first plot we show the frequencies in 
terms of the mass of the black hole for various values of Caoo . 
For each value of the speed of sound (each plot in Fig. [7|) we 
see that the bigger the velocity of the wind the bigger the 
frequency allowed. On the other hand, the three different 
panels show that the bigger the value of the speed of sound 
Csoo the bigger the frequency as well. 

In order to compare with observations, we include in 
Fig. [7] the masses and frequencies corresponding to the case 
of Sgr A*, where the mass of the supermassive black hole is 
M = 4.1±0.6x lO'^Mo (|Ghez et al. 2008|) and the frequency 
range was the one reported in (|Aschenbach et al. 20041 
lAbramowicz et al. 2004|) . We also show a shadowed region 
corresponding to the range of masses and frequencies asso- 
ciated to high mass X-ray binaries (HMXBs). 

Based on our numerical results in Fig. [7] we notice that 
for the range of velocities analyzed, for the model with 
r = 4/3 and various wind velocities, the frequencies of the 
density oscillations and mass of the black hole for Sgr A*, 
lie partially within that of observations. The trend observed 
between panels in Fig. [7| indicates that when Csoo increases, 
a bigger range of our parameters contain Sgr A*. 



Fig. [7] also shows a shadowed box corresponding to 
QPOs between ImHz to AOOmHz, observed in the spec- 
tra of HMXBs. Some of these measurements are associated 
with regions where it is believed a black hole exists, e.g Cyg 
X-1 (|Lui et al. 2004p . As we can see, this figure predicts a 
variety of possible configurations for the observed frequen- 
cies. For example, the allowed mass of the black hole that 
is hosted in the observed regions can run from hundreds to 
millions of solar masses for all the models presented here. 

The vertical dashed line in Fig.[7]corresponds to another 
QPO source of 1.27 Hz (|Kaur et al. 2007|l . Consistent black 
hole masses would be of about lO'^ to 10^ Mq. 

An exhaustive further study of the parameter space, 
would include information on the equation of state, addi- 
tional conditions involving non-adiabtic and radiative pro- 
cesses like in (|Zanotti et al. 20111 [Roedig et al. 20T2| in 
more realistic models. However it would be interesting as 
a first step the use of realistic wind velocities in simulations, 
which may involve the use of fish-eyed type of coordinates 
to describe the background space-time or the use of com- 
pactified coordinates with special slicing conditions on the 
space-time that penetrate the event horizon and approach 
infinity at once, that would allow us to study a bigger phys- 
ical domain in terms of the accretion radius, using modest 
computational resources, and associating bigger resolution 
in the regions near the black hole, as done in other already 
studied cases of the propagation of scalar waves onto black 
hole space-times where the future null infinity can be con- 
tained in the numerical domain l|Cruz-Osorio et al. 201ip . 



6 CONCLUSIONS AND DISCUSSION 

We have studied numerically the axisymmetric relativistic 
Bondi-Hoyle accretion of a supersonic ideal gas onto a fixed 
Schwarzschild background spacetime described using hori- 
zon penetrating coordinates. Our results reproduce the mor- 
phology of the distribution of the rest mass density of the 
gas obtained in the Newtonian regime and in previous rela- 
tivistic studies and verify that a stable shock cone forms. 

With our approach we measured the consistency in the 
rest mass density within the shock cone, which is an indi- 
cation that the excision implemented inside the black hole 
works fine with accuracy and convergence. Our treatment 
contributes to the solution of the Bondi-Hoyle accretion 
problem at the scale of the accretor scale length. This will 
provide a better understanding of the gas dynamics near the 
black hole. The restriction of this local type of study in as- 
trophysical scenarios at the moment is that one needs the 
use of a small accretion radius in order to cover a sufficiently 
large numerical domain, which in turn implies that the ve- 
locities that can be studied are higher than those observed 
or estimated by theoretical models. 

One of the properties of the accretion of a relativis- 
tic supersonic wind is that the density within the shock 
cone vibrates. We explore the possibility that such vibra- 
tions show frequencies within the range of those observed 
in QPO sources. Our approach is the first one for axisym- 
metric flows, in the relativistic regime, on a space-time cor- 
responding to a black hole, which additionally is described 
in coordinates that truly allow the gas flow inside the black 
hole horizon. 



Axisymmetric Bondi-Hoyle accretion onto a Schwarzschild Black Hole: shock cone vibrations 9 



10 I — 
1e-08 



HMXBs 
SgrA 

■XTE JOT 11. 2-731 7 

Mach=l 
Mach=2 
Mach=3 
Mach=4 



1e-06 0.0001 0.01 
F(Hz) 




HMXBs L--. 

SgrA m 

'XTE Ml 11. 2-731 7 ■ ■ 

Mach=1 
Mach=2 

Mach=3 - 

Mach=4 ■■■■ 



0.01 
F(Hz) 



HMXBs 
SgrA 

\XTEJ0U1.2-73n 

M3ch=1 
Mach=2 
Mach=3 
Mach=4 



0.0001 0.01 
F(Hz) 



Figure 7. We summarize our study of the parameter space. From 
top to bottom, the speed of sound used is = 0.1, 0.08, 0.05. 
We also show a vertical line indicating a special frequency of 
1.27Hz associated to HMXBs. We also indicate a shadowed region 
corresponding to HMXBs in the range of ImHz - 400mHz. The 
black box corresponds to the black hole mass and frequency range 
observed in Sgr A*. We show the case for F = 4/3. 



7 APPENDIX: TESTS OF THE CODE 

In order to illustrate how our implementation handles the 
evolution of a gas in a relativistic scenario, in this appendix 
we present standard tests. In all of them we use the minmod 
to reconstruct the variables and the HLLE flux formula, as 
we do in the simulations presented in this paper. 

Test 1, Id tests. The basic test a hydrodynamical code 
has to satisfy Sod's shock tube problem. We choose the ini- 
tial parameters of the commonly called blast wave problem. 



top- left 


top-right 


p = 0.1 


p = 0.1 


p = 1.0 


p = 0.01 


= 0.99 


= 


vy = 


yy = 


bottom-left 


bottom-right 


p = 0.5 


p = 0.1 


p = 1.0 


p = 1.0 


= 


v"^ =0 


vy = 


vy = 0.99 



Table 2. Initial data for the 2D shock. 



with standard parameters as in (|Marti' fc Miiller 2003 |l . The 
most stringent test has to do with the density profile at 
the front shock and the velocity of the shock. In Fig. [8] we 
show our results and contrast them with the exact solution, 
which we calculate following (|Martf fc Miiller 1994)) . The 
initial conditions are pL = pa = 1.0, pl = 1000, pr = 0.01, 

VL = Vr — 0. 

In order to test our implementation under a more ex- 
treme condition, we present what we call a head-on stream 
collision, which produces Lorentz factors of the order of 
1000. The initial conditions on the initial state are: pL ~ 
PR = 0.0009999998749, pi = pr = 3.33333258 x 10~^ 
vl = —Vr = 0.9999995, which are consistent with the wall- 
shock test in ( [Siegler fc Riffert 1999[ ). The numerical results 
compared with the exact solution are shown in Fig. [9] 

Test 2, 2d shocks. Our code handles the topology of 
cartesian coordinates because we use excision, which allows 
us to show standard 2d tests. The initial conditions are con- 
structed in such a way that a square domain is divided into 
four equal regions and where the variables start with the 
values shown in Tabled A snapshot at t = 0.35 is shown in 
Fig. 1101 We show the morphology as done in standard tests 
Pel Zanna fc Bucciantini 20021 IDel Zanna et al. 2007] . 

An additional test we want to include is the relativis- 
tic Emery's flat faced step. It consists in the flow entering 
from the left side of the domain and faces a step. The ini- 
tial conditions are p = 1.4, p = 1.0, = 0.99, u** = 
with F = 1.4. Three snapshots are shown in Fig. 1111 The 
boundary conditions are inflow at the left boundary, out- 
flow at the right, reflection at the top and bottom and the 
step boundaries. Again we show only the morphology as in 
Ponat et al. 1998p . 

Test 3, Michel accretion onto a black hole. The 
standard test on curved space-times is the radial ac- 
cretion of an ideal gas onto a Schwazschild black hole 
(|Michel 1972|) . For this, we constructed the exact solution in 
Eddington-Finkelstein coordinates using the parameters in 
( [Papadopoulos fc Font 1998] ) at coordinate time t — lOOM. 
The results are shown in Fig. 1121 where the distribution of 
the gas variables remains time-independent as expected. We 
also include a plot of the L2 norm of the error in time for 
various resolutions, showing the convergence of the applica- 
tion. 
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Figure 8. We show the gas variables at time t = 0.6 and the exact solution. As a global convergence test involves regions where constant 
states are involved, the velocity of the shock, and the shell of density at the front shock are important quantities to monitor. We resolve 
the shock with only a few points and the velocity at the time of the snapshot is very accurate. We use for this test 1000 cells and F = 5/3. 




Figure 12. Numerical and exact solution of the Michel accretion at time t = 0.990Af. The solid line is the exact solution and the dots 
are the numerical solution. The last plot shows the convergence order calculated using 100,200,400 cells; this shows the nearly second 
order convergence of our implementation of this test. The time is in units of the black hole mass M. 
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Figure 10. Logarithm of the rest mass density for the 2D rela- 
tivistic Riemann problem at t = 0.35. We illustrate the morphol- 
ogy developed by the process. We use 1200x1200 cells to cover 
the domain. 
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